"""
Phase 7: Structural Break Analysis
Tests whether demographic effects on rates and inflation changed around the GFC,
and whether QE and ZLB regimes moderate these effects.

Outputs:
  - Table 12: Pre-GFC vs post-GFC with Chow test (phase7_table12_chow.md)
  - Table 13: QE interaction (phase7_table13_qe.md)
  - Table 13b: QE + Japanification controls (phase7_table13b_japan_controls.md)
  - Table 13c: Japanification vs QE horse race (phase7_table13c_japan_horserace.md)
  - Table 14: ZLB interaction (phase7_table14_zlb.md)
  - Table 15: Rolling 15-year windows (phase7_table15_rolling.md)
  - Figure: Rolling Z_1 coefficient (output/figures/phase7_rolling_Z1_coef.png)
  - Data: Rolling Z_1 raw data (phase7_rolling_Z1.csv)
"""

import sys
from pathlib import Path
import numpy as np
import pandas as pd

# ── paths ────────────────────────────────────────────────────────────────────
PROJECT = Path("/mnt/c/demographics_capital_flows")
sys.path.insert(0, str(PROJECT / "multilateral" / "src"))
from model import PanelGLS

DATA_PATH = PROJECT / "monetary" / "data" / "processed" / "monetary_panel.csv"
TABLE_DIR = PROJECT / "monetary" / "output" / "tables"
FIG_DIR = PROJECT / "monetary" / "output" / "figures"
TABLE_DIR.mkdir(parents=True, exist_ok=True)
FIG_DIR.mkdir(parents=True, exist_ok=True)

OECD_38 = [
    "AUS", "AUT", "BEL", "CAN", "CHL", "COL", "CRI", "CZE", "DNK", "EST",
    "FIN", "FRA", "DEU", "GRC", "HUN", "ISL", "IRL", "ISR", "ITA", "JPN",
    "KOR", "LVA", "LTU", "LUX", "MEX", "NLD", "NZL", "NOR", "POL", "PRT",
    "SVK", "SVN", "ESP", "SWE", "CHE", "TUR", "GBR", "USA",
]


# ── helpers ──────────────────────────────────────────────────────────────────
def run_model(panel, dep_var, rhs_vars):
    """Run PanelGLS; return dict with coefs/se/pvals/r2/nobs or None."""
    cols = [dep_var] + rhs_vars + ["iso3", "year"]
    df = panel[cols].dropna()
    if len(df) < 50:
        print(f"  SKIP {dep_var}: only {len(df)} obs after dropna")
        return None
    try:
        y = df[dep_var].values
        X = df[rhs_vars].values
        gls = PanelGLS()
        gls.fit(y, X, df["iso3"].values, df["year"].values)
        resid = y - X @ gls.beta
        return {
            "coefs": dict(zip(rhs_vars, gls.beta)),
            "se": dict(zip(rhs_vars, gls.se)),
            "pvals": dict(zip(rhs_vars, gls.pvalues)),
            "r2": gls.r_squared,
            "nobs": gls.n_obs,
            "ncountries": gls.n_countries,
            "ssr": float(np.sum(resid ** 2)),
            "k": len(rhs_vars),
        }
    except Exception as e:
        print(f"  ERROR {dep_var}: {e}")
        return None


def star(p):
    """Significance stars."""
    if p < 0.01:
        return "***"
    elif p < 0.05:
        return "**"
    elif p < 0.10:
        return "*"
    return ""


def fmt_coef(res, var):
    """Format coefficient as 'coef(se)***'."""
    if res is None or var not in res["coefs"]:
        return ""
    c = res["coefs"][var]
    s = res["se"][var]
    p = res["pvals"][var]
    return f"{c:.3f}{star(p)} ({s:.3f})"


def fmt_val(val, fmt=".3f"):
    """Format a numeric value."""
    if val is None or np.isnan(val):
        return ""
    return f"{val:{fmt}}"


def build_markdown_table(title, row_vars, col_labels, results, footer_rows):
    """Build a pipe-format markdown table."""
    lines = [f"# {title}\n"]
    header = "| Variable | " + " | ".join(col_labels) + " |"
    sep = "|" + "---|" * (len(col_labels) + 1)
    lines.append(header)
    lines.append(sep)
    for var in row_vars:
        cells = [fmt_coef(r, var) for r in results]
        lines.append(f"| {var} | " + " | ".join(cells) + " |")
    lines.append(sep)
    for label, vals in footer_rows:
        lines.append(f"| {label} | " + " | ".join(vals) + " |")
    return "\n".join(lines)


def chow_test(ssr_pooled, ssr_1, ssr_2, n1, n2, k):
    """Compute Chow test F-statistic and p-value.

    F = ((SSR_pooled - SSR_1 - SSR_2) / k) / ((SSR_1 + SSR_2) / (n1 + n2 - 2k))
    """
    from scipy import stats as sp_stats

    df1 = k
    df2 = n1 + n2 - 2 * k
    if df2 <= 0:
        return np.nan, np.nan
    F = ((ssr_pooled - ssr_1 - ssr_2) / df1) / ((ssr_1 + ssr_2) / df2)
    p = 1 - sp_stats.f.cdf(F, df1, df2)
    return F, p


# ── main ─────────────────────────────────────────────────────────────────────
def main():
    print("=" * 70)
    print("PHASE 7: STRUCTURAL BREAK ANALYSIS")
    print("=" * 70)

    panel = pd.read_csv(DATA_PATH)
    print(f"Panel: {len(panel):,} obs, {panel['iso3'].nunique()} countries")

    oecd = panel[panel["iso3"].isin(OECD_38)].copy()
    pre_gfc = panel[panel["year"] <= 2007].copy()
    post_gfc = panel[panel["year"] >= 2008].copy()

    print(f"OECD:     {len(oecd):,} obs, {oecd['iso3'].nunique()} countries")
    print(f"Pre-GFC:  {len(pre_gfc):,} obs ({pre_gfc['year'].min()}-{pre_gfc['year'].max()})")
    print(f"Post-GFC: {len(post_gfc):,} obs ({post_gfc['year'].min()}-{post_gfc['year'].max()})")

    z_vars = ["Z_1", "Z_2", "Z_3"]
    controls_rate = ["rgdp_growth", "inflation", "fiscal_bal_gdp", "kaopen", "nfa_gdp_lag"]
    controls_infl = ["rgdp_growth", "output_gap", "fiscal_bal_gdp", "kaopen", "nfa_gdp_lag"]

    # ── TABLE 12: Pre-GFC vs Post-GFC with Chow Test ─────────────────────
    print("\n── Table 12: Pre-GFC vs Post-GFC (Chow Test) ──")

    dvs_12 = [
        ("real_bond_10y", controls_rate),
        ("inflation", controls_infl),
    ]

    results_t12 = []
    col_labels_t12 = []
    chow_rows = []

    for dv, controls in dvs_12:
        rhs = z_vars + controls

        # Pooled
        res_pool = run_model(panel, dv, rhs)
        # Pre-GFC
        res_pre = run_model(pre_gfc, dv, rhs)
        results_t12.append(res_pre)
        col_labels_t12.append(f"{dv} (Pre-GFC)")
        # Post-GFC
        res_post = run_model(post_gfc, dv, rhs)
        results_t12.append(res_post)
        col_labels_t12.append(f"{dv} (Post-GFC)")

        # Chow test
        if res_pool and res_pre and res_post:
            F, p = chow_test(
                res_pool["ssr"], res_pre["ssr"], res_post["ssr"],
                res_pre["nobs"], res_post["nobs"], res_pool["k"],
            )
            chow_rows.append((dv, F, p))
            print(f"  {dv}: Chow F = {F:.3f}, p = {p:.4f}")
            print(f"    Pre-GFC Z_1 = {res_pre['coefs']['Z_1']:.3f} "
                  f"(p={res_pre['pvals']['Z_1']:.3f})")
            print(f"    Post-GFC Z_1 = {res_post['coefs']['Z_1']:.3f} "
                  f"(p={res_post['pvals']['Z_1']:.3f})")
        else:
            chow_rows.append((dv, np.nan, np.nan))

    row_vars_t12 = list(dict.fromkeys(z_vars + controls_rate + controls_infl))

    footer_t12 = [
        ("R²", [f"{r['r2']:.3f}" if r else "" for r in results_t12]),
        ("N obs", [f"{r['nobs']:,}" if r else "" for r in results_t12]),
        ("N countries", [f"{r['ncountries']}" if r else "" for r in results_t12]),
    ]
    # Add Chow statistics as footer rows
    for dv, F, p in chow_rows:
        # Put Chow stat under the two columns belonging to this DV
        chow_vals = [f"F={F:.3f}" if not np.isnan(F) else "", f"p={p:.4f}" if not np.isnan(p) else ""]
        # Pad to full width
        idx = next(i for i, (d, _, _) in enumerate(chow_rows) if d == dv) * 2
        vals = [""] * len(results_t12)
        vals[idx] = chow_vals[0]
        vals[idx + 1] = chow_vals[1]
        footer_t12.append((f"Chow ({dv})", vals))

    md_t12 = build_markdown_table(
        "Table 12: Pre-GFC vs Post-GFC Structural Break",
        row_vars_t12, col_labels_t12, results_t12, footer_t12,
    )
    print("\n" + md_t12)
    t12_path = TABLE_DIR / "phase7_table12_chow.md"
    t12_path.write_text(md_t12)
    print(f"\nSaved: {t12_path}")

    # ── TABLE 13: QE Interaction ──────────────────────────────────────────
    print("\n── Table 13: QE Interaction ──")

    # Ensure interaction terms exist
    for z in z_vars:
        col = f"{z}_x_qe"
        if col not in panel.columns:
            panel[col] = panel[z] * panel["qe_active"]

    qe_interactions = ["Z_1_x_qe", "Z_2_x_qe", "Z_3_x_qe"]
    qe_countries = panel[panel["qe_country"] == 1].copy()
    non_qe_countries = panel[panel["qe_country"] == 0].copy()

    # Ensure interaction terms in subsets
    for z in z_vars:
        col = f"{z}_x_qe"
        if col not in qe_countries.columns:
            qe_countries[col] = qe_countries[z] * qe_countries["qe_active"]
        if col not in non_qe_countries.columns:
            non_qe_countries[col] = non_qe_countries[z] * non_qe_countries["qe_active"]

    print(f"QE countries:     {qe_countries['iso3'].nunique()} countries, "
          f"{len(qe_countries):,} obs")
    print(f"Non-QE countries: {non_qe_countries['iso3'].nunique()} countries, "
          f"{len(non_qe_countries):,} obs")

    dvs_13 = [
        ("real_bond_10y", controls_rate),
        ("inflation", controls_infl),
    ]

    results_t13 = []
    col_labels_t13 = []

    for dv, controls in dvs_13:
        # Full sample with QE interactions
        rhs_int = z_vars + ["qe_active"] + qe_interactions + controls
        res = run_model(panel, dv, rhs_int)
        results_t13.append(res)
        col_labels_t13.append(f"{dv} (Full+QE int)")
        if res:
            for iv in qe_interactions:
                print(f"  {dv} Full: {iv}={res['coefs'][iv]:.3f} "
                      f"(p={res['pvals'][iv]:.3f})")

        # QE countries only
        rhs_qe = z_vars + ["qe_active"] + controls
        res = run_model(qe_countries, dv, rhs_qe)
        results_t13.append(res)
        col_labels_t13.append(f"{dv} (QE only)")
        if res:
            print(f"  {dv} QE-only: Z_1={res['coefs']['Z_1']:.3f} "
                  f"(p={res['pvals']['Z_1']:.3f})")

        # Non-QE countries only
        rhs_noqe = z_vars + controls
        res = run_model(non_qe_countries, dv, rhs_noqe)
        results_t13.append(res)
        col_labels_t13.append(f"{dv} (Non-QE)")
        if res:
            print(f"  {dv} Non-QE: Z_1={res['coefs']['Z_1']:.3f} "
                  f"(p={res['pvals']['Z_1']:.3f})")

    row_vars_t13 = list(dict.fromkeys(
        z_vars + ["qe_active"] + qe_interactions + controls_rate + controls_infl
    ))

    footer_t13 = [
        ("R²", [f"{r['r2']:.3f}" if r else "" for r in results_t13]),
        ("N obs", [f"{r['nobs']:,}" if r else "" for r in results_t13]),
        ("N countries", [f"{r['ncountries']}" if r else "" for r in results_t13]),
    ]

    md_t13 = build_markdown_table(
        "Table 13: QE Interaction with Demographics",
        row_vars_t13, col_labels_t13, results_t13, footer_t13,
    )
    print("\n" + md_t13)
    t13_path = TABLE_DIR / "phase7_table13_qe.md"
    t13_path.write_text(md_t13)
    print(f"\nSaved: {t13_path}")

    # ── TABLE 14: ZLB Interaction ─────────────────────────────────────────
    print("\n── Table 14: ZLB Interaction ──")

    zlb_var = "Z_1_x_near_zlb"
    if zlb_var not in panel.columns:
        panel[zlb_var] = panel["Z_1"] * panel["near_zlb"]

    dvs_14 = [
        ("real_bond_10y", controls_rate),
        ("inflation", controls_infl),
    ]

    results_t14 = []
    col_labels_t14 = []

    for dv, controls in dvs_14:
        # Full sample with ZLB interaction
        rhs_zlb = z_vars + ["near_zlb", zlb_var] + controls
        res = run_model(panel, dv, rhs_zlb)
        results_t14.append(res)
        col_labels_t14.append(f"{dv} (Full)")
        if res and zlb_var in res["coefs"]:
            print(f"  {dv} Full: Z_1_x_near_zlb={res['coefs'][zlb_var]:.3f} "
                  f"(p={res['pvals'][zlb_var]:.3f})")

        # OECD
        # Ensure interaction in oecd
        if zlb_var not in oecd.columns:
            oecd[zlb_var] = oecd["Z_1"] * oecd["near_zlb"]
        res = run_model(oecd, dv, rhs_zlb)
        results_t14.append(res)
        col_labels_t14.append(f"{dv} (OECD)")
        if res and zlb_var in res["coefs"]:
            print(f"  {dv} OECD: Z_1_x_near_zlb={res['coefs'][zlb_var]:.3f} "
                  f"(p={res['pvals'][zlb_var]:.3f})")

    row_vars_t14 = list(dict.fromkeys(
        z_vars + ["near_zlb", zlb_var] + controls_rate + controls_infl
    ))

    footer_t14 = [
        ("R²", [f"{r['r2']:.3f}" if r else "" for r in results_t14]),
        ("N obs", [f"{r['nobs']:,}" if r else "" for r in results_t14]),
        ("N countries", [f"{r['ncountries']}" if r else "" for r in results_t14]),
    ]

    md_t14 = build_markdown_table(
        "Table 14: ZLB Interaction with Demographics",
        row_vars_t14, col_labels_t14, results_t14, footer_t14,
    )
    print("\n" + md_t14)
    t14_path = TABLE_DIR / "phase7_table14_zlb.md"
    t14_path.write_text(md_t14)
    print(f"\nSaved: {t14_path}")

    # ── TABLE 15: Rolling 15-Year Windows ─────────────────────────────────
    print("\n── Table 15: Rolling 15-Year Windows ──")

    rolling_records = []

    for start_year in range(1985, 2011):
        end_year = start_year + 14
        window = panel[(panel["year"] >= start_year) & (panel["year"] <= end_year)].copy()
        mid_year = start_year + 7  # label for the window

        for dv, controls, dv_label in [
            ("real_bond_10y", controls_rate, "10y_rate"),
            ("inflation", controls_infl, "inflation"),
        ]:
            rhs = z_vars + controls
            res = run_model(window, dv, rhs)
            if res:
                rolling_records.append({
                    "start_year": start_year,
                    "end_year": end_year,
                    "mid_year": mid_year,
                    "dv": dv_label,
                    "Z_1_coef": res["coefs"]["Z_1"],
                    "Z_1_se": res["se"]["Z_1"],
                    "Z_1_pval": res["pvals"]["Z_1"],
                    "Z_2_coef": res["coefs"]["Z_2"],
                    "Z_2_se": res["se"]["Z_2"],
                    "Z_2_pval": res["pvals"]["Z_2"],
                    "Z_3_coef": res["coefs"]["Z_3"],
                    "Z_3_se": res["se"]["Z_3"],
                    "Z_3_pval": res["pvals"]["Z_3"],
                    "r2": res["r2"],
                    "nobs": res["nobs"],
                    "ncountries": res["ncountries"],
                })

    rolling_df = pd.DataFrame(rolling_records)
    csv_path = TABLE_DIR / "phase7_rolling_Z1.csv"
    rolling_df.to_csv(csv_path, index=False)
    print(f"Saved rolling data: {csv_path}")
    print(f"  {len(rolling_df)} window-DV combinations")

    # Build summary markdown table
    md_lines = ["# Table 15: Rolling 15-Year Window Z_1 Coefficients\n"]
    md_lines.append("| Window | DV | Z_1 coef | Z_1 SE | Z_1 p-val | R² | N |")
    md_lines.append("|--------|----|---------:|-------:|----------:|---:|--:|")

    for _, row in rolling_df.iterrows():
        s = star(row["Z_1_pval"])
        md_lines.append(
            f"| {int(row['start_year'])}-{int(row['end_year'])} "
            f"| {row['dv']} "
            f"| {row['Z_1_coef']:.3f}{s} "
            f"| {row['Z_1_se']:.3f} "
            f"| {row['Z_1_pval']:.4f} "
            f"| {row['r2']:.3f} "
            f"| {int(row['nobs']):,} |"
        )

    md_t15 = "\n".join(md_lines)
    t15_path = TABLE_DIR / "phase7_table15_rolling.md"
    t15_path.write_text(md_t15)
    print(f"Saved: {t15_path}")

    # ── Figure: Rolling Z_1 coefficient with confidence bands ─────────
    print("\n── Figure: Rolling Z_1 Coefficient ──")
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

    for ax, dv_label, title in [
        (axes[0], "10y_rate", "Real 10-Year Bond Yield"),
        (axes[1], "inflation", "Inflation"),
    ]:
        sub = rolling_df[rolling_df["dv"] == dv_label].sort_values("mid_year")
        if len(sub) == 0:
            continue

        mid = sub["mid_year"].values
        coef = sub["Z_1_coef"].values
        se = sub["Z_1_se"].values
        upper = coef + 1.96 * se
        lower = coef - 1.96 * se

        ax.plot(mid, coef, "b-o", markersize=4, linewidth=1.5, label="Z_1 coefficient")
        ax.fill_between(mid, lower, upper, alpha=0.2, color="blue", label="95% CI")
        ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")

        # Shade QE period (2008-2022) — show as shaded band on midpoint axis
        qe_start = max(mid.min(), 2008 - 7)  # midpoint when 2008 enters window
        qe_end = min(mid.max(), 2022 - 7)    # midpoint when 2022 exits window
        # Actually shade 2008-era midpoints: windows centered around GFC
        ax.axvspan(2004, 2017, alpha=0.1, color="red", label="QE era windows")

        ax.set_ylabel("Z_1 Coefficient")
        ax.set_title(f"Rolling 15-Year Z_1 on {title}")
        ax.legend(fontsize=8, loc="best")
        ax.grid(True, alpha=0.3)

    axes[1].set_xlabel("Window Midpoint Year")
    plt.tight_layout()

    fig_path = FIG_DIR / "phase7_rolling_Z1_coef.png"
    plt.savefig(fig_path, dpi=150, bbox_inches="tight")
    plt.close()
    print(f"Saved: {fig_path}")

    # ── TABLE 13b: Japanification Controls in Structural Break ─────────
    print("\n── Table 13b: Japanification Controls in Structural Break ──")

    # Horse race: does Z₁×QE survive after controlling for Japanification proxies?
    # Proxies: rgdp_growth (expected growth proxy), debt_liab_gdp, savings_investment_gap
    japan_controls_rate = [
        "rgdp_growth", "inflation", "fiscal_bal_gdp", "kaopen", "nfa_gdp_lag",
        "debt_liab_gdp", "savings_investment_gap",
    ]

    # Ensure QE interactions exist in panel
    for z in z_vars:
        col = f"{z}_x_qe"
        if col not in panel.columns:
            panel[col] = panel[z] * panel["qe_active"]

    qe_interactions = ["Z_1_x_qe", "Z_2_x_qe", "Z_3_x_qe"]

    results_t13b = []
    col_labels_t13b = []

    # Col 1: Full sample — QE interaction + Japanification controls
    rhs_13b_full = z_vars + ["qe_active"] + qe_interactions + japan_controls_rate
    res = run_model(panel, "real_bond_10y", rhs_13b_full)
    results_t13b.append(res)
    col_labels_t13b.append("10y (Full+Japan)")
    if res:
        print(f"  Full+Japan: Z_1_x_qe={res['coefs']['Z_1_x_qe']:.3f} "
              f"(p={res['pvals']['Z_1_x_qe']:.3f})")
        print(f"    Z_1={res['coefs']['Z_1']:.3f} (p={res['pvals']['Z_1']:.3f})")

    # Col 2: OECD — QE interaction + Japanification controls
    for z in z_vars:
        col = f"{z}_x_qe"
        if col not in oecd.columns:
            oecd[col] = oecd[z] * oecd["qe_active"]
    res = run_model(oecd, "real_bond_10y", rhs_13b_full)
    results_t13b.append(res)
    col_labels_t13b.append("10y (OECD+Japan)")
    if res:
        print(f"  OECD+Japan: Z_1_x_qe={res['coefs']['Z_1_x_qe']:.3f} "
              f"(p={res['pvals']['Z_1_x_qe']:.3f})")

    # Col 3: Full sample — QE interaction WITHOUT Japanification (baseline comparison)
    rhs_13b_base = z_vars + ["qe_active"] + qe_interactions + controls_rate
    res = run_model(panel, "real_bond_10y", rhs_13b_base)
    results_t13b.append(res)
    col_labels_t13b.append("10y (Full, no Japan)")

    # Col 4: OECD — QE interaction WITHOUT Japanification (baseline comparison)
    res = run_model(oecd, "real_bond_10y", rhs_13b_base)
    results_t13b.append(res)
    col_labels_t13b.append("10y (OECD, no Japan)")

    row_vars_t13b = list(dict.fromkeys(
        z_vars + ["qe_active"] + qe_interactions + japan_controls_rate
    ))

    footer_t13b = [
        ("R²", [f"{r['r2']:.3f}" if r else "" for r in results_t13b]),
        ("N obs", [f"{r['nobs']:,}" if r else "" for r in results_t13b]),
        ("N countries", [f"{r['ncountries']}" if r else "" for r in results_t13b]),
        ("Japanification controls", ["Yes", "Yes", "No", "No"]),
    ]

    md_t13b = build_markdown_table(
        "Table 13b: QE Interaction with Japanification Controls",
        row_vars_t13b, col_labels_t13b, results_t13b, footer_t13b,
    )
    print("\n" + md_t13b)
    t13b_path = TABLE_DIR / "phase7_table13b_japan_controls.md"
    t13b_path.write_text(md_t13b)
    print(f"\nSaved: {t13b_path}")

    # ── TABLE 13c: Direct Japanification Interactions ─────────────────
    print("\n── Table 13c: Direct Japanification Interactions ──")

    # Create binary Japanification indicators
    median_growth = panel["rgdp_growth"].median()
    median_debt = panel["debt_liab_gdp"].median()
    panel["low_growth"] = (panel["rgdp_growth"] < median_growth).astype(float)
    panel["high_debt"] = (panel["debt_liab_gdp"] > median_debt).astype(float)
    panel["Z_1_x_low_growth"] = panel["Z_1"] * panel["low_growth"]
    panel["Z_1_x_high_debt"] = panel["Z_1"] * panel["high_debt"]

    oecd["low_growth"] = (oecd["rgdp_growth"] < median_growth).astype(float)
    oecd["high_debt"] = (oecd["debt_liab_gdp"] > median_debt).astype(float)
    oecd["Z_1_x_low_growth"] = oecd["Z_1"] * oecd["low_growth"]
    oecd["Z_1_x_high_debt"] = oecd["Z_1"] * oecd["high_debt"]

    japan_int_vars = ["low_growth", "high_debt", "Z_1_x_low_growth", "Z_1_x_high_debt"]

    results_t13c = []
    col_labels_t13c = []

    # Col 1: Full — Japanification interactions only
    rhs_13c_japan = z_vars + japan_int_vars + controls_rate
    res = run_model(panel, "real_bond_10y", rhs_13c_japan)
    results_t13c.append(res)
    col_labels_t13c.append("10y (Full, Japan int)")
    if res:
        print(f"  Full Japan-only: Z_1_x_low_growth={res['coefs']['Z_1_x_low_growth']:.3f} "
              f"(p={res['pvals']['Z_1_x_low_growth']:.3f})")
        print(f"    Z_1_x_high_debt={res['coefs']['Z_1_x_high_debt']:.3f} "
              f"(p={res['pvals']['Z_1_x_high_debt']:.3f})")

    # Col 2: OECD — Japanification interactions only
    res = run_model(oecd, "real_bond_10y", rhs_13c_japan)
    results_t13c.append(res)
    col_labels_t13c.append("10y (OECD, Japan int)")
    if res:
        print(f"  OECD Japan-only: Z_1_x_low_growth={res['coefs']['Z_1_x_low_growth']:.3f} "
              f"(p={res['pvals']['Z_1_x_low_growth']:.3f})")

    # Col 3: Full — Both QE and Japanification interactions (horse race)
    rhs_13c_both = z_vars + ["qe_active"] + qe_interactions + japan_int_vars + controls_rate
    res = run_model(panel, "real_bond_10y", rhs_13c_both)
    results_t13c.append(res)
    col_labels_t13c.append("10y (Full, both)")
    if res:
        print(f"  Full both: Z_1_x_qe={res['coefs']['Z_1_x_qe']:.3f} "
              f"(p={res['pvals']['Z_1_x_qe']:.3f})")
        print(f"    Z_1_x_low_growth={res['coefs']['Z_1_x_low_growth']:.3f} "
              f"(p={res['pvals']['Z_1_x_low_growth']:.3f})")

    # Col 4: OECD — Both QE and Japanification interactions (horse race)
    res = run_model(oecd, "real_bond_10y", rhs_13c_both)
    results_t13c.append(res)
    col_labels_t13c.append("10y (OECD, both)")
    if res:
        print(f"  OECD both: Z_1_x_qe={res['coefs']['Z_1_x_qe']:.3f} "
              f"(p={res['pvals']['Z_1_x_qe']:.3f})")

    row_vars_t13c = list(dict.fromkeys(
        z_vars + ["qe_active"] + qe_interactions + japan_int_vars + controls_rate
    ))

    footer_t13c = [
        ("R²", [f"{r['r2']:.3f}" if r else "" for r in results_t13c]),
        ("N obs", [f"{r['nobs']:,}" if r else "" for r in results_t13c]),
        ("N countries", [f"{r['ncountries']}" if r else "" for r in results_t13c]),
        ("QE interactions", ["No", "No", "Yes", "Yes"]),
        ("Japanification interactions", ["Yes", "Yes", "Yes", "Yes"]),
    ]

    md_t13c = build_markdown_table(
        "Table 13c: Japanification vs QE — Horse Race",
        row_vars_t13c, col_labels_t13c, results_t13c, footer_t13c,
    )
    print("\n" + md_t13c)
    t13c_path = TABLE_DIR / "phase7_table13c_japan_horserace.md"
    t13c_path.write_text(md_t13c)
    print(f"\nSaved: {t13c_path}")

    # ── Summary ──────────────────────────────────────────────────────────
    print("\n── Key Findings ──")

    # Chow test results
    for dv, F, p in chow_rows:
        sig = "SIGNIFICANT structural break" if p < 0.05 else "no significant break"
        print(f"  Chow test ({dv}): F={F:.3f}, p={p:.4f} -> {sig}")

    # QE moderation
    res_qe_rate = results_t13[0]
    if res_qe_rate and "Z_1_x_qe" in res_qe_rate["coefs"]:
        c = res_qe_rate["coefs"]["Z_1_x_qe"]
        p = res_qe_rate["pvals"]["Z_1_x_qe"]
        print(f"  QE × Z_1 on rates: {c:.3f} (p={p:.3f})")

    # ZLB moderation
    res_zlb_rate = results_t14[0]
    if res_zlb_rate and zlb_var in res_zlb_rate["coefs"]:
        c = res_zlb_rate["coefs"][zlb_var]
        p = res_zlb_rate["pvals"][zlb_var]
        print(f"  ZLB × Z_1 on rates: {c:.3f} (p={p:.3f})")

    # Rolling stability
    rate_rolling = rolling_df[rolling_df["dv"] == "10y_rate"]
    if len(rate_rolling) > 0:
        sig_pct = (rate_rolling["Z_1_pval"] < 0.10).mean() * 100
        print(f"  Rolling (10y rate): Z_1 significant in {sig_pct:.0f}% of windows")

    print("\nPhase 7 complete.")


if __name__ == "__main__":
    main()
